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1.  Introduction 


Under  the  sponsorship  of  the  Defense  Nuclear  Agency  (DNA),  the  Physics 
Internationa  Company  developed  a  Transportable  Electromagnetic  Pulse 
Simulator  (TEMPS)  for  Harry  Diamond  Laboratories  (HDL)*  in  1972.  In 
1974,  HDL  procured  a  copy  of  TEMPS,  called  AESOP  (Army  electromag¬ 
netic  pulse  (EMP)  Simulator  Operation),  which  was  erected  at  the 
Woodbridge  Research  Facility  (WRF)  (see  fig.  1).  In  1994,  AESOP  was 
transferred  to  the  Test  and  Evaluation  Command  (TECOM)  at  White 
Sands  Missile  Range.  The  main  function  of  AESOP  is  to  simulate  the  elec¬ 
tromagnetic  envirorunent  produced  by  an  exoatmospheric  nuclear  burst. 
By  simulating  the  threat-level  EMP  environment,  it  allows  for  testing  and 
verifying  the  EMP  susceptibility  or  hardness  of  a  system.  Since  AESOP  is  a 
virtual  copy  of  TEMPS  (except  for  transportability  and  some  minor  modi¬ 
fications  and  enhancements),  the  specifications  of  AESOP  are  documented 
in  the  TEMPS's  technical  manuals. 

Because  of  the  complexity  of  EMP  simulation  and  the  limited  computing 
capability  in  the  70s,  no  computing  codes  were  developed  to  calculate  or 
approximate  the  field  generated  by  AESOP.  A  field  mapping  report^  of 
TEMPS  was  developed  to  provide  the  EMP  fields  at  key  locations  near  the 
ground.  The  present  report  provides  results  from  a  code  that  calculates  the 
electromagnetic  (EM)  field  generated  by  AESOP  (the  code  is  given  in 
app  A). 


Figure  1.  AESOP  at  WRF. 


*Noiu  part  of  the  Army  Research  Laboratory. 

^Eugene  L.  Patrick  and  Spencer  L.  Soohoo,  Transportable  Electromagnetic  Pulse  Simulator  (TEMPS)  Field  Map¬ 
ping  Report,  Harry  Diamond  Laboratories,  HDL-TR-1743  (Reprinted  May  1984). 


2.  AESOP  Simulator 


Most  of  the  information  in  this  section  is  derived  from  the  Physics  Interna¬ 
tional  manual  for  TEMPS.^ 

Located  at  WRF  in  a  120,000-m^  test  area,  AESOP  is  the  largest  EMP  free- 
field  simulator  in  the  world.  It  consists  of  three  basic  subsystems;  a  sup¬ 
port  structure,  pulse  generator,  and  antenna. 

•  The  support  structure,  made  of  fiberglass  columns  and  laminated  wood 
beams,  is  used  to  support  and  suspend  the  pulse  generator  and  antenna  at 
a  height  up  to  20  m  above  the  ground. 

•  The  pulse  generator  mainly  comprises  the  pulser  assembly  with  two  Marx 
generators,  peaking  capacitors,  output  switch,  and  control  subsystem. 

•  The  antenna  is  composed  of  a  cylindrical  and  conical  system  of  36  wires 
extending  from  each  end  of  the  pulser  assembly  with  associated  forming 
hoops,  end  assemblies,  tensioning  elements,  and  end  grounding  termina¬ 
tions.  The  diameter  of  the  cylindrical  section  is  9.144  m. 

The  general  layout  of  the  AESOP  system  is  shown  in  figure  2.^  The  overall 
length  of  AESOP  is  300  m  from  end  dome  to  end  dome.  The  pulser  and 
antenna  are  normally  positioned  at  20  m  above  the  ground  level  but  can 
also  be  at  lower  elevations. 


^Physics  International  Company,  Transportable  Electromagnetic  Pulse  Simulator  System  (TEMPS)  Operation  & 
Maintermnce  Manual,  PIMM-372  (April  1973). 
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AESOP  produces  a  pulsed  EM  field  by  transmitting  a  high-energy  electri¬ 
cal  pulse  generated  from  the  pulse  generator  to  the  groimd  through  the 
anteima.  It  can  provide  a  radiated  free-field  range  of  20  to  50  kV/m  from 
the  pulser  on  the  center  line.  The  pulse  shape  is  a  double  exponential  with 
less  than  10-percent  undershoot  (fig.  3).  For  more  information,  see  the 
TEMPS  manual.^ 


Figure  3.  AESOP 
characteristics; 

(a)  pulse  shape  and 

(b)  span  on  EMP 
coverage. 


3.  Theory 


In  the  model  presented  here,  the  AESOP  anterma  is  modeled  as  a  series  of 
dipole  antenna  elements  above  a  perfectly  conducting  groimd.  There  are 
two  steps  involved  in  the  calculation  of  the  radiated  electromagnetic  field 
from  AESOP.  A  closed-form  analytical  time-domain  current  profile  on 
each  dipole  element  is  derived  by  the  transmission-line  method.  (In  order 
to  simplify  the  problem,  we  accoimt  only  for  reflection  at  the  load  end  of 
AESOP.)  The  time-domain  radiated  electromagnetic  field  from  each  dipole 
is  then  computed  by  a  convolution  technique. 


^Physics  International  Company,  Transportable  Electromagnetic  Pulse  Simulator  System  (TEMPS)  Operation  & 
Maintenance  Manual,  PIMM-372  (April  1973). 
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3.1  Current  Distribution 

AESOP  is  represented  by  a  transmission  line  over  a  perfectly  conducting 
ground.  There  are  three  factors  in  this  transmission-line  model  that  need  to 
be  considered.^ 

First,  the  35-stage  Marx  generator  can  be  modeled  in  a  lumped  equivalent 
circuit,  as  shown  in  figure  4(a).  The  switches  and  Sj  are  the  voltage 
charging  and  discharging  switches,  respectively.  is  the  characteristic 
impedance  of  the  antenna,  which  includes  the  biconical  and  conical  sec¬ 
tions.  Zj  is  the  termination  load. 

Second,  the  biconical  section  with  cones  at  angle  6  =  40.5°  has  a  character¬ 
istic  impedance,  Zj, , 


=  [cot  (I)]  =  120  £2  , 


(1) 


Figure  4.  Circuit 
models:  (a)  AESOP, 

(b)  transmission-line 
model,  and  (c)  AESOP 
model  at  t  =  ^2* 
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^Physics  International  Company,  Transportable  Electromagnetic  Pulse  Simulator  System  (TEMPS)  Operation  & 
Maintenance  Manual,  PIMM-372  (April  1973). 
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where 


• 

Third,  the  cylindrical  system  of  wires  that  makes  up  the  AESOP  antenna 
structure  can  be  treated  as  low-loss  transmission  line  of  a  wire  over  the 
ground  plane.  Its  characteristic  impedance,  Zc,  is  approximated  by  R,  L, 
and  C  (R  is  series  resistance  per  unit  length,  L  is  series  inductance  per  unit 
length,  and  C  is  shunt  capacitance  per  unit  length): 


where 


L 


s  =  jco, 

(0  =  angular  frequency, 

h  =  height  of  a  wire  over  the  groimd  plane,  and 
a  =  radius  of  a  wire. 

The  conical  section  is  designed  to  maintain  the  same  characteristic  imped¬ 
ance  as  the  cylindrical  section.  In  other  words,  the  ratio  of  2/i/fl  remains 
constant  for  both  conical  and  cylindrical  sections. 

For  the  low-resistive-loss  biconical  section,  the  characteristic  impedance, 
Zb,  is  very  compatible  with  the  characteristic  impedance,  Zc,  of  die  cylin¬ 
drical  section.  Since  the  cylindrical  section  is  much  longer  than  the 
biconical  section,  the  assumption  is  that  the  characteristic  impedance  Zb 
equals  the  characteristic  impedance  Zc.  The  schematic  representation  of  a 
1-cm-unit-length  section,  6i,  of  Zc  is  shown  in  figiure  4(b). 

At  time  t  =  t2,  we  obtain  the  Laplace  transform  for  the  AESOP  model  (fig. 
4(c)): 


Vo  =  Vit2) 


[  7  \ 

Wc 

[ft, (-»,)] 

(3) 
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where 


^1/  ^2/  ^3  are  initial  voltage  sources  at  t  =  / 

Zeq  is  the  equivalent  impedance  of  R2,  L2,  Zb,  Zq,  and  Z/, 
yi/  VI  are  roots  of 


'^(^2)  Ri 

1/  -U 

^1(^2)  ,  1  1 

m2)C2  £-1 

y  + 

^(^2)  ^1^2 

1/  -A 

^3  2L  ' 

^1/  52/  53, 54,  S5  are  roots  of 

s5  +  +  A2s3  +  A3S2  +  A4S  +  A5  =  0  ,  and 

Ai,  A-i,  A3,  A4,  As  are  functions  of  R,  R^,  R2,  L,  Li,  L2,  C,  Ci,  C2. 


The  reflecting  voltage  coefficient  at  the  load  end  is 


Voltage  across  the  capacitor,  V„,  at  any  given  section  of  the  antenna  is 
calculated  as 

^n  =  ^ok”^+r'e-(2^^-«)V^  ,  (5) 

where 


yf=a  +  jP  , 


The  corresponding  current  on  the  «th  section  of  the  antenna  is 


I  L) 


(6) 


The  time-domain  current,  /„(f),  can  be  calculated  directly  from  the  above 
equation  by  the  Laplace  transform  technique: 


I„{t)  =  [B„it  exp  {D„^t)]  u(t  -  Q  ,  (7) 

where 

Bfi](  is  the  time  constant  coefficient  number  for  the  A:th  variable  and  the 
section, 

Dfik  is  the  time  constant  exponential  coefficient  number  for  the  variable 
and  the  section,  and 


tn  is  the  delay  time  constant  for  the  section. 

The  calculated  current  distributions  at  various  locations  on  the  AESOP  an¬ 
tenna  are  shown  in  figure  5(a).  Figure  5(b)  compares  the  results  of  this 
method  to  those  of  the  finite  time-difference  method.  Good  agreement  was 
obtained  at  the  early  time,  and  reasonable  agreement  was  obtained  at  the 
later  time. 


3.2  Electric  Field  and  Magnetic  Field  of  a  Dipole 


Figtu-e  5.  Current 
comparisons:  (a)  on 
different  locations 
away  from  center  of 
AESOP  and  (b)  on 
bicone  and  at  1.3  ft 
away  from  center. 


As  described  in  section  2,  an  electromagnetic  field  is  emitted  from  the 
AESOP  antenna,  which  is  a  transmission-line  antenna.  It  can  be  modeled 
as  a  series  of  linear  dipoles.  Each  dipole  itself  is  an  elementary  anteima  that 
radiates  fields.  In  this  section,  the  electric  field  and  magnetic  field  (radiated 
free  field)  equations  of  a  dipole  model  are  expanded  from  section  2.3  of 
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In  the  more  general  case  shown  in  figure  6(b),  r{x,y;z)  and  r'{x',y' ,z')  are  the 
coordinates  of  the  observation  point  and  of  the  dipole,  respectively.  Hence, 
R{X,Y,Z)  is  the  relative  coordinate  of  the  observation  point  in  the  Cartesian 
coordinate  with  the  origin  at  the  dipole,  where  R=r  -7'.  The  magnitude 
of  R  can  be  defined  as 

R  =  \/(x-xf4.(!/-yf+(z-zf  . 

In  the  general  case,  then,  equations  (8)  and  (9)  become 


E(o))  =  ^(X  +  -^ 


In 


^(oR- 


ZQdU  j(o 


An 


cR 


I{a))e-l^^^(=cosdR 

— I{o})  sin  0  0  , 

jcoR^ 


(10) 


H{co)  =  ^  ^  +  “V]  sin  0  O  . 

ycK 


(11) 


By  transferring  equations  (10)  and  (11)  into  the  time  domain  and  changing 
to  Cartesian  coordinates,  we  transform  equations  (10)  and  (11)  to 


An 


R 


2R" 


XZ 

R^ 


Zodl 


47t 


2R'^ 


^9 


3Z^-R^ 


R^ 


J 


(12) 


z  , 


R  R^J  ' 


(13) 


where  sgn  W  “  j|j  '  ®  is  a  convolution  operator. 

We  can  change  the  orientation  of  the  dipole  so  that  it  matches  those  of 
AESOP  by  rotating  the  coordinates  as  indicated  in  figure  6(c); 

from  x  to  y 
y  to  2 
2  to  AC 
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Finally,  the  electric  and  magnetic  fields  of  a  dipole  are  calculated  by 


PI  A  -  ZlO 


E(f)  = 


Zndl 


Idllf  R]X^-R^  . 
cRdt\  c}  ^2 


R 


3X^-R\ 


R^ 


Zodl 


iTT 

Zodl 


iTT 


Mg 


xz  = 

R^ 


(14) 


cRdt 


+ 


X 

r} 


(15) 


Implementation 

The  dipole  current  distribution  and  field  equations  developed  in  section  3 
provide  us  the  means  to  compute  the  free  field  emitted  by  AESOP.  But  the 
total  field  in  the  half  space  above  the  grovmd  is  a  result  of  the  superposi¬ 
tion  of  the  free  field  and  the  reflected  field  from  the  ground. 

To  simplify  the  problem,  we  assume  that  the  earth  grotmd  is  a  perfectly 
electrically  .conducting  (PEC)  half  space.  With  this  assumption,  weean  ap¬ 
ply  the  method  of  images  to  calculate  the  reflected  field.  The  calculation 
can  be  further  simplified  if  we  use  the  symmetry  of  dipole  currents 
through  the  plane  that  is  perpendicular  to  the  axis  of  AESOP  at  the  center. 
The  biconic  and  cylindrical  parts  of  the  antenna  are  assumed  to  be  one 
single  transmission  line,  having  both  ends  slanted  toward  the  ground. 

We  model  the  simulator  as  a  series  of  dipole  elements  representing  the 
physical  structure  of  the  AESOP  simulator.  The  number  of  dipole  elements 
has  been  selected  to  resolve  the  early-time  structure  of  the  field  and  to  pro¬ 
duce  a  relatively  smooth  waveform.  For  every  0.01-m  section  of  the  an¬ 
tenna,  we  evaluate  the  current  distribution  developed  in  section  3.1  at  the 
time  and  spatial  location  that  is  the  driving  source  for  the  corresponding 
dipole  element.  We  employ  14,978  dipole  elements  with  a  dipole  separa¬ 
tion  of  0.01  m  to  represent  half  the  antenna  of  AESOP. 

We  let  the  coordinates  be  as  depicted  in  figure  7.  Thus,  the  z  =  0  plane  is 
the  ground  plane  and  parallel  to  AESOP.  The  method  of  images  is  applied 
through  the  ground  plane  to  compute  the  reflected  fields.  The  y  =  0  plane 
contains  the  antenna.  The  symmetry  of  current  distribution  is  through  the 
X  =  0  plane,  called  the  center  plane.  The  height  of  the  antenna  is  20  m  above 
the  ground.  Each  end  of  AESOP  is  terminated  at  6.5  m  above  the  groimd. 
Each  slanted  section  begins  and  ends  at  100.4  and  147.9  m  from  the  center 
plane,  respectively.  To  save  computational  time,  we  translate  and  rotate 
the  coordinate  system  to  determine  the  fields  of  the  simulator  slanted  sec- 


tions,  as  depicted  in  figure  8.  Finally,  we  use  superposition  of  all  fields  to 
obtain  the  total  field  at  an  observation  point  from  the  AESOP  simulator. 

We  have  used  a  FORTRAN  program  to  compute  all  fields  at  an  observa¬ 
tion  point.  The  program  is  listed  in  appendix  A. 


z 
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5.  Results  and  Comparison  with  Empirical  Data 

We  implemented  a  FORTRAN  program  to  compute  the  fields  at  an  obser¬ 
vation  point  out  to  2  ps  in  time  with  a  time  step  of  1  ns.  The  following  are 
the  execution  times  of  the  program  on  some  systems: 

System  Execution  time  (hourrminute) 

IBM  3090  with  MVS/XA  2.3  4:00 

IBM  RISC/6000-530  with  AIX  3.2.4  3:00 

IBM  RISC/6000-550  with  AIX  3.2.4  2:18 

IBM  RISC/6000-560  with  AIX  3.2.4  1:50 

SGI  4D/440  with  IRIX  4.0.5  2:32 

SGI  4D/35  with  IRIX  4.0.5  11:11 

DECstation  5000  with  ULTRIX  4.3  14:49 

Because  of  the  symmetry,  the  computed  field  components  Ex,  Hy,  and  Hz 
are  the  only  nonzero  components  of  the  field  on  the  center  plane  (x  =  0 
plane).  Off  the  center  plane,  the  computed  model  indicates  that  the  Hx 
component  is  zero  (in  retarded  time)  imtil  the  current  arrives  to  the  slanted 
section  of  the  simulator.  In  the  plane  containing  the  antenna  (y  =  0  plane), 
the  symmetry  requires  that  Ex,  Ez,  and  Hy  be  the  only  nonzero  fields.  The 
calculated  fields  show  the  effects  of  the  change  in  direction  of  the  current 
on  the  antenna  structure,  as  slope  discontinuities  in  the  fields,  at  the 
expected  time. 

For  comparison  to  the  empirical  data,  we  scale  all  the  computed  fields  by  a 
factor  that  scales  the  computed  total  electric  field  of  the  point  (0,50,20)  to 
50  kV /m.  In  comparison,  we  notice  that  the  idealized  lossless  ground 
plane  assumption  has  caused  the  peak  fields  of  the  model  to  be  generally 
overestimated  (compared  to  the  empirical  data)  by  roughly  a  factor  of  two; 
this  assumption  also  led  to  the  inaccuracy  of  the  vertical  field  prediction 
from  the  slanted  end  sections  of  the  simulator.  Assuming  a  perfectly  con¬ 
ducting  grovmd  plane  means  that  the  model  does  not  reflect  the  effects  of 
the  pseudo-Brewster  angle  on  the  vertical  fields.  However,  the  energy  in 
the  predicted  field  waveforms  still  provides  an  upper  bound  for  the  actual 
fields. 

Appendix  B  contains  field  plots  compared  to  empirical  data  for  the  follow¬ 
ing  sample  points  in  the  space: 
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1546, 

4.6) 

(410, 

502, 

4.6) 

(-444, 

1341, 

4.6) 

(  473, 

0, 

1) 

(  473, 

0, 

4.6) 
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6.  Conclusion 


Although  the  perfectly  electrically  conducting  ground  assumption  has 
generally  overestimated  and  distorted  the  calculated  field,  the  predicted 
fields  are  still  good  for  bounding  problems,  especially  for  points  high 
above  the  groimd.  The  calculated  information  can  provide  upper  boimd 
contours  or  isosurfaces  at  certain  field  magnitudes.  The  number  of  dipoles 
used  in  the  model  can  be  decreased,  with  a  corresponding  decrease  in  the 
model  run  time,  but  with  an  increase  in  "spiky"  noise  in  the  predicted 
waveforms.  This  model  has  been  helpful  in  imderstanding  the  relationship 
between  the  characteristics  of  the  field  waveforms  and  the  structure  of  the 
AESOP  simulator.  Because  it  is  simple  and  requires  only  2.5  MB  of 
memory  to  execute,  this  model  can  be  used  on  most  computer  platforms. 
This  code  is  the  first  building  block  to  develop  more  accurate  field  predic¬ 
tions  for  the  AESOP  simulator. 
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Appendix  A. — Program  Listing 


IMPLICIT  NONE 

INTEGER*4  MAXF,  MAXS,  MAXT,  MAXSH 

PARAMETER  (MAXF=8,  MAXS=14978,  MAXSH=10040,  MAXT=2000) 
CHARACTER  DSCUCO*80,  DSEHD*80,  DSEHT*80,  SPC*1 
REAL*8  COEF(20,MAXS) ,  CUC 
REAL* 8  CUCC,  CUL,  CUR 

REAL*8  EH(0:MAXT,MAXF) ,  RL(4),  RO(12),  RS(4) 

REAL*8  EHM(8),  R0S(4),  ROT{4) 

REAL*8  OXT,  OZT,  SA,  SAX,  SAZ,  SINA,  COSA 
REAL*8  AL,  C,  DL,  DLMIN,  DT,  PI,  ZO,  Z004PI 
REAL*8  CU,  CUD,  CVO,  ETC,  HTC,  TCS,  TEHS,  TR,  RETARD 
REAL*8  CUP,  TCRS,  SN,  SNl,  SNR,  SNRl 

INTEGER*4  IB,  INF,  IT,  ITS,  IS,  ISB,  ISE,  ISL,  ISGX,  ISGZ 
INTEGER*4  I,  J,  K,  ITM(8),  ISGXB,  ISGXE,  ISGZB,  ISGZE 
COMMON  /COEFBL/COEF,CUC 

COMMON  /CUPARM/CUP, TCRS , SN, SNl , SNR, SNRl 

DATA  CUCC,CUL,CUR/0.273448D-12,4.07D-9,5.364D-3/ 

* 

CUP  =  DSQRT( CUL* CUCC) 

CUC  =  CUR/ (2 . *DSQRT(CUL/CUCC) ) 

* 

*  ■  INITIALIZE  PARAMETERS 

* 

ISGXB  =  1 
ISGXE  =  -1 
ISGZB  =  1 
ISGZE  =  -1 
SPC  = 

C  C  =  2.997925D8 

C  =  3.D8 

PI  =  4.*DATAN(1.D0) 

ZO  =  4.*PI*C*l.D-7 
C  ZO  =  376.991184D0 

Z004PI=C*l.D-7 
DLMIN  =  l.D-2 
DT  =  l.D-9 
AL  =  20. 

RS(2)  =  0. 

OXT  =  100.40D0 
OZT  =  AL 
SAX  =  47.5D0 
SAZ  =  13.5D0 

SA  =  DSQRT(SAX**2  +  SAZ* *2) 


COSA  =  SAX/SA 
SINA  =  SAZ/SA 
ISE  =  0 


READ  IN  CURRENT  COEFFICIENTS  OF  ALL  SEGMENTS 

DSCUCO  =  ' /home/emsc3/bluu/aesop/CUCOEF.DATA' 

OPEN (28,  FILE=DSCUCO ,  STATUS= ' OLD ' ,  FORM= ' FORMATTED ' ) 

DO  170  J=l,  MAXS 

READ(28,*)  IS,  (COEF{I,J),  1=1,20) 

170  CONTINUE 
CLOSE (28) 

READ  PARTITION  DATA  SET  NAME  OF  THE  E  AND  H  FIELD 

READ(5, ' (A) ' )  DSEHD 
IB  =  INDEX (DSEHD, '  ' ) 

IB  =  IB  -  1 

READ  IN  SEGMENT  LENGTH  AND  COORDINATE  OF  OBSERVATION  POINT 

115  READ(5, *,END=999)  DL,  (RO(I),  1=1,12) 

ISL  =  (DL.  +  DLMIN/2) /DLMIN 

DETERMINE  THE  STARTING  TIME 

ISB  =  (ISL+1) /2 
ISGX  =  1 

IF  (RO(l) .LT.O. )  ISGX  =  -1 
RS(1)  =  ISGX* (2 *ISB-1) * (DLMIN/2. ) 

TEHS  =  DSQRT( (RO(l) -RS (1) ) **2  +  RO(2)**2  +  (RO (3 ) -AL) **2 ) /C 
I  ISB*CUP 

INITIALIZE  EH  ARRAY 

DO  150  J=l,  MAXF 
DO  150  1=0,  MAXT 
EH(I,J)  =  0. 

150  CONTINUE 

CALCULATE  E  AND  H  FIELD 

WRITE (6,*)  'P',  RO(l),  RO(2),  RO(3) 

WRITE (6,*)  '  DIPOLE  LENGTH  IS',  DL 


DO  200  IS=ISB,  MAXSH,  ISL 
ISE  =  IS 
TCS  =  IS*CUP 
TORS  =  (2* (MAXS-IS)+1) *CUP 
SN  =  IS-1 
SNl  =  IS 

SNR  =  2*MAXS  -  IS 
SNR1=  2*MAXS  -  IS  +1 
DO  220  ISGZ=  ISGZB, ISGZE, -2 
DO  220  ISGX=  ISGXB, ISGXE, -2 

RS(1)  =  ISGX*(2*IS-l)*(DLMIN/2.) 

RS{3)  =  ISGZ*AL 

C  RS(4)  =  DSQRT(RS(1) **2  +  RS(2)**2  +  RS(3)**2) 

RL(1)  =  ROd)  -  RS(1) 

RL(2)  =  RO{2)  -  RS(2) 

RL(3)  =  RO{3)  -  RS(3) 

RL(4)  =  DSQRT(RL{1) **2  +  RL(2)**2  +  RL(3)**2) 

* 

RETARD  =  RL(4)/C  +  TCS  -  TEHS 
ITS  =  RETARD/DT  +  1. 

DO  250  IT=ITS,  MAXT 
TR  =  IT*DT  -  RETARD 
CALL  CURRENT (IS, TR,CU, CUD, C VO) 

HTC  =  ISGZ*DL^ (CUD/RL(4) /C  + 

I  CU/RL(4) **2) / (4.*PI) 

*  EH(IT,4)  =  EH(IT,4)  +  0. 

EH(IT,5)  =  EH(IT,5)  -  RL (3 ) /RL (4 ) *HTC 
EH(IT,6)  =  EH(IT,6)  +  RL (2 ) /RL (4 ) *HTC 
ETC  =  ISGZ*Z004PI*DL* (CUD/RL(4) /C  + 

I  3.*{CU/RL(4)**2)  +  (1.5*0  * (CVO/RL(4) **3) ) 

EH(IT,1)  =  EH(IT,1)  +  ISGZ*Z004PI*DL* ( 

I  CUD/RL(4) /C* ( (RL(1) /RL(4) ) **2-1. )  + 

I  (CU/RL(4)**2+(C/2.)*(CVO/RL(4)**3) )* 

I  (3.*(RL(1)/RL(4))**2-1.)  ) 

EH(IT,2)  =  EH(IT,2)  +  RL ( 1 ) /RL (4 ) *RL (2 ) /RL ( 4 ) *ETC 
EH(IT,3)  =  EH(IT,3)  +  RL ( 1 ) /RL (4 ) *RL (3 ) /RL ( 4 ) *ETC 
250  CONTINUE 

220  CONTINUE 
200  CONTINUE 

* 

*  CALCULATE  THE  FIELD  CONTRIBUTE  BY  THE  SLANT  PARTS 

* 

RS(3)  =  0. 

ROT (2)  =  RO(2) 

ISB  =  ISE  +  ISL 
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DO  260  IS=ISB,  MAXS,  ISL 
TCS  =  IS*CUP 

TORS  =  (2* (MAXS-IS) +1) *CUP 
SN  =  IS-1 
SNl  =  IS 

SNR  =  2*MAXS  -  IS 
SNR1=  2*MAXS  -  IS  +1 
DO  280  ISGZ=  ISGZB, ISGZE, -2 
DO  280  ISGX=  ISGXB, ISGXE, -2 
ROS(l)  =  R0{1)  -  ISGX*OXT 
ROS(3)  =  RO(3)  -  ISGZ*OZT 

R0T{1)  =  COSA*ROS(l)  -  ISGX*ISGZ*SINA*ROS (3 ) 

ROT(3)  =  ISGX*ISGZ*SINA*ROS (1)  +  COSA*ROS(3) 

RS(1)  =  ISGX* (2* (IS-MAXSH) -1) * {DLMIN/2. ) 

C  RS(4)  =  DSQRT(RS (1) **2  +  RS(2)**2  +  RS(3)**2) 

RL(1)  =  ROT(l)  -  RS{1) 

RL(2)  =  ROT (2)  -  RS(2) 

RL(3)  =  ROT(3)  -  RS(3) 

RL(4)  =  DSQRT(RL(1) **2  +  RL(2)**2  +  RL(3)**2) 

* 

RETARD  =  RL(4)/C  +  TCS  -  TEHS 
ITS  =  RETARD/DT  +  1. 

DO  290  IT=ITS,  MAXT 
TR  =  IT*DT  -  RETARD 
CALL  CURRENT ( I S , TR , CU , CUD , CVO ) 

ir 

HTC  =  ISGZ*DL* (CUD/RL(4) /C  + 

I  CU/RL(4) **2) / (4. *PI) 

EHM(6)  =  RL(2) /RL(4) *HTC 

★ 

EH(IT,4)  =  EH(IT,4)  +  ISGX*ISGZ*SINA*EHM { 6 ) 

EH(IT,5)  =  EH(IT,5)  -  RL ( 3 ) /RL (4 ) *HTC 
EH(IT,6)  =  EH(IT,6)  +  COSA*EHM(6) 

★ 

ETC  =  ISGZ*Z004PI*DL* (CUD/RL(4) /C  + 

I  3.*(CU/RL(4) **2)  +  (1.5*0  * (CV0/RL(4) **3) ) 

EHM(l)  =  ISGZ*Z004PI*DL* ( 

I  CUD/RL(4)/C*((RL(1)/RL(4))**2-1.)  + 

I  (CU/RL(4) **2+ (C/2 . ) * (CVO/RL(4) **3) ) * 

I  (3.*(RL(1)/RL(4) )**2-l.)  ) 

EHM(3)  =  RL(1) /RL(4) *RL(3) /RL(4) *ETC 

EH(IT,1)  =  EH(IT,1)  +  COSA*EHM(l) +ISGX*ISGZ*SINA*EHM(3) 
EH(IT,2)  =  EH(IT,2)  +  RL ( 1 ) /RL (4 ) *RL ( 2 ) /RL (4 ) *ETC 
EH(IT,3)  =  EH(IT,3)  -  ISGX*ISGZ*SINA*EHM ( 1 ) +C0SA*EP1M (3 ) 
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290  CONTINUE 

280  CONTINUE 
260  CONTINUE 


DO  270  1=1,8 
EHM(I)  =  0. 

ITM{I)  =  0. 

270  CONTINUE 

DETERMINE  THE  MAX  VALUES  OF  E  &  H  FIELDS  AND  THEIR  TIMES 
DO  295  IT=0,  MAXT 

ETC  =  DSQRT{EH(IT, 1) **2  +  EH(IT,2)**2  +  EH(IT,3)**2) 
HTC  =  DSQRT(EH(IT, 4) **2  +  EH(IT,5)**2  +  EH(IT,6)**2) 
EH(IT,7)=ETC 
EH(IT, 8)=HTC 
IF  (ETC.GT.EHM{7) )  THEN 
EHM(7)  =  ETC 
ITM(7)  =  IT 
ENDIF 

IF  (HTC.GT.EHM(8) )  THEN 
EHM(8)  =  HTC 
.  ITM(8)  =  IT 
ENDIF 

DO  297  1=1,6 

IF  (ABS(EH(IT,I) ) .GT.ABS(EHM{I) ) )  THEN 
EHM{I)  =  EH (IT, I) 

ITM(I)  =  IT 
ENDIF 

297  CONTINUE 
295  CONTINUE 

SET  A  DATA  SET  NAME  FOR  E  AND  H  FIELD 

WRITE (DSEHT, ’ (SP, A, E14 . 7 , E14 . 7 , E14 . 7 ) ' )  DSEHD(1:IB) , 

I  RO(l),  RO(2),  RO(3) 

IF ( DSEHT ( IB+1 : IB+1 ) .NE. •  ’)  THEN 

DSEHT (IB+2:IB+2)  =  DSEHT (IB+1 : IB+1 ) 

DSEHT (IB+16:IB+16)  =  DSEHT ( IB+15 : IB+15 ) 
DSEHT(IB+30:IB+30)  =  DSEHT(IB+29 : IB+29) 

ENDIF 

DSEHT (IB+1: IB+1)  =  SPC 

DSEHT (IB+15: IB+15)  =  SPC 
DSEHT (IB+29: IB+29)  =  SPC 


RETARD  =  TEHS/DT 


WRITE  E  AND  H  FIELD  TO  FILE 

OPEN (32,  FILE=DSEHT,  STATUS= ' UNKNOWN ' ,  FORM= ' UNFORMATTED ' ) 
WRITE(32)  RO(l) , RO (2 ) , RO (3 ) , RO (4 ) ,RO(5) ,RO(6) ,DL,DT 
WRITE ( 32 )  RO ( 7 ) , RO ( 8 ) , RO ( 9 ) , RO ( 10 ) , RO ( 11 ) , RO ( 12 ) , TEHS , RETARD 
WRITE (32)  EHM(l) ,EHM(2) ,EHM(3) ,EHM(4) ,EHM(5) ,EHM(6) ,EHM(7) , 

I  EHM(8) 

WRITE(32)  ITM(l) ,ITM(2) ,ITM(3) ,ITM(4) ,ITM(5) ,ITM(6) , 

I  ITM(7) , ITM(8) 

DO  400  1=0,  MAXT 

WRITE(32)  (EH(I,J),  J=l,  8) 

400  CONTINUE 
CLOSE (32) 

WRITE (6,*)  'END  OF  THIS  RUN;  DATA  IN  FILE:  ',  DSEHT 
GOTO  115 
999  STOP 

79  FORMAT(8E15.7) 

END 


SUBROUTINE  CALCULATE,  THE  DIPOLE ' S 

CURRENT, 

CURRENT  DERIVATIVE, 

CONVOLUTION  OF  SIGN  FUNCTION  AND  CURRENT 


SUBROUTINE  CURRENT  ( NS ,  TS ,  Y ,  YD ,  YCO ) 

IMPLICIT  NONE 

INTEGER *4  MAXS ,  NS,  LR 

REAL*8  P,  TC2,  SN,  SNl ,  SNR,  SNRl 

REAL*8  T,  TS,  Y,  YS,  YD,  YSD,  YCO,  YSCO 

COMMON  /CUPARM/P, TC2 , SN, SNl , SNR, SNRl 

DATA  P/0.333606558688524581D-10/ 

DATA  C/0.219836066437990879D-04/ 


LR  =  0 
T  =  TS 

CALL  CURDIS (NS, LR, SNl , T, YS , YSD, YSCO) 
T  =  TS  +  P 

CALL  CURDIS (NS, LR,SN,T,Y, YD, YCO) 

Y  =  Y  -  YS 


YD  =  YD  -  YSD 
YCO  =  YCO  -  YSCO 

* 

*  REFLECTION  CURRENT  IS  NOT  EMERGING 

* 

IF  (T.LE.TC2)  RETURN 
T  =  TS  -  TC2 

* 

*  REFLECTION  CURRENT 

* 

LR  =  12 

CALL  CURDIS (NS , LR, SNRl , T , YS , YSD, YSCO) 

Y  =  Y  +  YS 
YD  =  YD  +  YSD 
YCO  =  YCO  +  YSCO 
T  =  T  +  P 

CALL  CURDIS (NS , LR, SNR, T, YS , YSD, YSCO) 

Y  =  Y  -  YS 
YD  =  YD  -  YSD 
YCO  =  YCO  -  YSCO 
RETURN 

END 

* 

* 

**  SUBROUTINE  CALCULATES  CURRENT  PARAMETERS  ** 
* 

★ 

SUBROUTINE  CURDIS (NS, IR, SN, T, YS, YSD, YSCO) 
IMPLICIT  NONE 
INTEGER *4  MAXS,  NS,  IR 
PARAMETER  (MAXS=14978) 

REAL*8  CO(20,MAXS),  C 
REAL* 8  T,  SN,  YS,  YSD,  YSCO 
COMMON  /COEFBL/CO,C 


* 

AA 

= 

CO( 

1+IR,NS) 

■k 

AB 

= 

CO{ 

2+IR,NS) 

* 

AC 

= 

CO( 

3+IR,NS) 

★ 

AD 

= 

CO( 

4+IR,NS) 

* 

PD 

= 

CO( 

5+IR,NS) 

* 

AF 

= 

CO{ 

6+IR,NS) 

* 

EA 

= 

CO( 

7, NS) 

★ 

EB 

= 

CO( 

8, NS) 

★ 

EC 

= 

CO( 

9. NS) 

★ 

ED 

= 

CO (10, NS) 
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★ 

★ 

* 

* 

* 

* 

* 

* 

* 

★ 

★ 

★ 

★ 

★ 


EDP 

EF 

RA 

RB 

RC 

RD 

PDR 

RF 

AR 

ER 


C0(11,NS) 
CO (12, NS) 
CO(13,NS) 
CO (14, NS) 
CO(15,NS) 
CO (16, NS) 
CO (17, NS) 
CO (18, NS) 
CO (19, NS) 
CO(20,NS) 


CURRENT  DISTRIBUTION 

YS  =  DEXP(-SN*C) * (  CO(l+IR,NS) *DEXP(CO(7,NS) *T)  + 

I  CO(2+IR,NS) *DEXP(CO(8,NS) *T)  + 

I  CO(3+IR,NS) *DEXP(CO(9,NS) *T)  + 

I  CO(6+IR,NS) *DEXP(CO(12,NS) *T)  + 

I  2 .*CO(4+IR,NS) *DEXP(CO(10,NS) *T) *DCOS (CO (11,NS) *T+CO(5+IR,NS) )  ) 

DERIVATIVE  OF  CURRENT  DISTRIBUTION 


YSD  =  DEXP(-SN*C) * (  CO ( 1+IR, NS) *CO ( 7 , NS) *DEXP (CO ( 7 , NS ) *T)  + 

CO(2+IR,NS) *CO(8,NS) *nEXP(CO(8,NS) *T)  + 
CO(3+IR,NS) *CO(9,NS) *DEXP (CO ( 9 , NS ) *T)  + 
CO(6+IR,NS) *CO(12,NS) *DEXP (CO ( 12 , NS) *T) 
2 . *CO(4+IR,NS) *DEXP(CO(10,NS) *T) * ( 

CO (10, NS) *DCOS(CO (11,NS) *T+CO ( 5+IR, NS ) ) - 
CO(ll,NS) *DSIN(CO(ll,NS) *T+CO (5+IR,NS) ) )  ) 


CONVOLUTION  OF  SGN(T)  AND  I(T) 


YSCO  =  DEXP(-SN*C) * ( 

CO (1+IR, NS) * (2 . *DEXP(CO(7,NS) *T) -1. ) /CO (7, NS)  + 
CO(2+IR,NS) * (2. *DEXP(CO(8,NS) *T) -1 . ) /CO ( 8 , NS)  + 
CO(3+IR,NS) * (2 . *DEXP(CO(9,NS) *T) -1. ) /CO (9, NS)  + 
CO(6+IR,NS) * (2. *DEXP(CO(12,NS) *T) -1 . ) /CO ( 12 , NS )  + 

2 . *CO(4+IR,NS) * (  2 . *DEXP(CO(10,NS) *T) * 

(CO(10,NS) *DCOS (CO(ll,NS) *T+CO ( 5+IR, NS ) )+ 
CO(ll,NS) *DSIN(CO(ll,NS) *T+CO ( 5+IR, NS ) ) ) 
(CO(10,NS)*DCOS(CO(5+IR,NS) )+ 

CO(ll,NS) *DSIN(CO(5+IR,NS) ) )  )/ 
(CO(10,NS) **2+CO(ll,NS) **2)  ) 

IF(IR.EQ.O)  RETURN 

YS=  YS  +  DEXP(-SN*C) *CO(19,NS) *DEXP(CO(20,NS) *T) 
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YSD=  YSD  +  DEXP(-SN*C) *C0(19,NS) *CO(20,NS) *DEXP(CO(20,NS) *T) 
YSCO=  YSCO  +  DEXP(-SN*C)*CO(19,NS)*{2.*DEXP(CO(20,NS)*T)-1.)/ 
CO{20,NS) 

RETURN 

END 
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Appendix  B. — Field  Plot  Comparison 

The  field  plots  presented  display  electric  or  magnetic  field  component 
waveforms  at  selected  test  points,  showing  both  empirical  measxirements 
and  computed  results  in  the  same  graph.  No  comparisons  are  done  for 
field  components  that  are  zero  in  theory.  Both  of  the  fields,  empirical  and 
computed,  are  scaled  to  the  same  reference  field  level,  50  kV/m  at  50  m 
away  from  the  center  of  AESOP  on  the  center  plane.  We  selected  the  field 
points  based  on  the  avahabihty  of  empirical  data  for  test  points  that  were 
measured  for  the  Perimeter  Field  Map  at  the  Woodbridge  Research  Facil¬ 
ity.  For  those  test  points  that  have  more  than  one  measurement,  we  have 
chosen  the  best  data  set  to  compare  with  computed  data.  The  selected  field 
comparison  points  are  as  follows: 
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Figure  B-1  (cont'd).  (c) 
Field  comparison  at 
(0,-176,1)  meters: 


Figure  B-3.  Field 
comparison  at 
(-332,"1764)  meters: 
(a)  E^,  (b)  Ey. 
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Figure  B-3  (confd). 
Field  comparison  at 
(-332,-176,1)  meters: 
(c)Ej,(d)H^. 
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(d) 


Figure  B-3  (conrd).  (e) 
Field  comparison  at 
(-332,-176,1)  meters: 

(e)Hy,(£)H^. 
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:  Data  scaled  to  50  kV/m 
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Figure  B-4.  Field  (a) 
comparison  at 
(-332,-176,4.6)  meters:  ^ 
(a)  Ex,  (b)  Ey. 
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Figure  B-4  (cont'd).  (c) 
Field  comparison  at 
(-332-176,4.6)  meters:  ^ 

(c)£2,(d)H,. 
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Figure  B-4  (cont'd).  (e) 
Field  comparison  at 
(-332,-176,4.6)  meters:  ® 


(e)Hy,(f)H2. 
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Figure  B-5  (cont'd). 
Field  comparison  at 
(0,1546,4-6)  meters: 
(c)H,. 
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Figure  B-6.  Field 
comparison  at 
(410,502,4*6)  meters: 
(a)E^,  (b)Ey. 
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Figure  B-6  (cont'd). 
Field  comparison  at 
(410,502,4.6)  meters: 
(c)E^,(d)H,. 
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Figure  B-6  (cont'd). 
Field  comparison  at 
(410,502,4-6)  meters: 
(e)Hy,(f)H^. 
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Figure  B-7.  Field  (a) 
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Figure  B-7  (cant'd).  (e) 
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Figure  B-8.  Field  (a) 
comparison  at  (473,0,1) 
meters:  (a)  E^,  (b)  E^.  | 
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Figure  B-8  (cont'd). 
Field  comparison  at 
(473A1)  meters:  (c)  Hy. 


(C) 
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Figure  B-9.  Field 
comparison  at 
(473,0,4.6)  meters: 
(a)  E^,  (b)  E^. 


^Empirical 


Computed 


Data  scaled  to  50  kV/m 

— i - i - 1 _ i _ i _ i _ i _ i _ i _ 

200  400  600  800  1000  1200  1400  1600  1800  2C 

Time  (ns) 


200  400  600  800  1000  1200  1400  1600  1800  2( 

Time  (ns) 


Figure  B-9  (cant'd).  (c) 

Field  comparison  at 

(473,0,4*6)  meters: 

(c)Hy. 
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